Testing
Passes all CTests on pm-cpu and pm-gpu, except for IO_TEST. I'm not sure why this test would be affected by this PR.

Passes Polaris ocean/horiz_press_grad/salinity_gradient, ocean/horiz_press_grad/temperature_gradient, and ocean/horiz_press_grad/ztilde_gradient. Convergence plots below:

ocean/horiz_press_grad/salinity_gradient

ocean/horiz_press_grad/temperature_gradient

ocean/horiz_press_grad/ztilde_gradient

Omega pressure gradient has 10^-16 - 10^-13 differences with the python implementation in polaris.

Polaris omega_pr suite

• Baseline workdir: /pscratch/sd/s/sbrus/polaris_pgrad_pr_suite_baseline
• Baseline build: /global/homes/s/sbrus/scratch/polaris_pgrad_pr_suite_baseline/build
• PR build: /global/homes/s/sbrus/scratch/polaris_pgrad_pr_suite/build
• PR workdir: /global/homes/s/sbrus/scratch/polaris_pgrad_pr_suite
• Machine: pm-cpu
• Compiler: gnu
• Build type: Release
• Log: not found
• Result:
  • Failures (2 of 9):
    • ocean/spherical/qu/rotation_2d
    • ocean/spherical/icos/cosine_bell/restart
  • Diffs (1 of 9):
    • ocean/horiz_press_grad/salinity_gradient

ocean/spherical/qu/rotation_2d also failed in baseline.
I'm looking into ocean/spherical/icos/cosine_bell/restart fail.

The IOSTREAM_TEST only fails on a re-run when the test/ocn.*.nc files are already present. I noticed this also happens in the develop branch.

The ocean/spherical/icos/cosine_bell/restart polaris test fail was fixed by @cbegeman's suggestion to turn off the pressure gradient tendency. I will submit a follow-on polaris PR that does that for cases that do not require a pressure gradient term.

The IOTest fail was an issue on my end with having an out-of-date scorpio submodule.

This now passes CTests with gnu Release/Debug on pm-cpu and gnugpu Release/Debug on pm-gpu.

Thanks for the feedback @mwarusz and @brian-oneill. I think I've addressed your comments in the the last couple commits.

@mwarusz, I choose to add TracerTimeLevel to the computeVelocityTendenciesOnly, computeVelocityTendencies, and computeAllTendencies. Let me know if you think that works okay or if you'd like to see changes.

Related to the TracerTimeLevel change is where we actually want to call computeSpecVol. I'm open to moving it to another location.

This passes CTests on pm-gpu and pm-cpu as well as the omega_pr suite on pm-cpu.

@sbrus89

Oops, turns out that using TracerTimeLevel is not really doing the right thing in RungeKutta4 currently. I'm sorry for suggesting that change. The reason why it is not correct is because these calls

Array2DReal Temp     = Tracers::getByName(TracerTimeLevel, "Temperature");
Array2DReal Salinity = Tracers::getByName(TracerTimeLevel, "Salinity");

get temperature and salinity from the tracers array in the singleton Tracers instance. However, intermediate RK4 stages should get these tracers from the ProvisTracers array.

In general, I believe that we should have TracerTimeLevel. However, in order to make it work, we would need to refactor the Tracers class to not be a singleton. Then ProvisTracers could become a proper Tracers class instance, instead of a plain array.

@mwarusz , that makes sense with respect to RK4. I'm open to making Tracers a non-singleton class, but I suppose we could also pass the 3D Tracer array (either from Tracers::getAll() or the ProvisTracers array) and use IndxTemp and IndxSalt to access the right subviews. But that essentially duplicates what the Tracers::getByIndex() function does, so maybe it's better to go the non-singleton route.

@mwarusz, I am now passing a tracer array into computeVelocityTendenciesOnly(). I left the TracerTimeLevel argument even though it isn't serving a functional purpose at this point. I can remove it or leave it for now until we switch Tracers to a non-singleton class.

Tested head on perlmutter gnu CPU, GPU, and passes all cTests.

I tested the baroclinic channel using polaris_0.10.0-alpha.1 and the setup

polaris setup -t  ocean/planar/baroclinic_channel/1km/default --model mpas-ocean -p $e/master/components/mpas-ocean -w $r/260401-BclCh-MpAB2

etc for both MPAS-O and Omega, for resolutions 10km, 4km, 1km. All use the default namelist except I set horizontal viscosity (del2) to 20 m^2/s^2 to match the Petersen et al. 2015 paper.

Omega simulations generally have the right behavior, but with some notable differences. Here is

• left: Omega 1km RK4, dt=2s nu=20
• right: MPAS. 1km AB2, dt=30s nu=20

https://github.com/user-attachments/assets/c19c2d18-381e-46cc-9b8e-a021dc6f1588
Video might not work. Here are snapshots:
Day 3: Something is wrong on the temperature gradients early on. Perhaps horizontal advection is anti-diffusive, or vertical advection is moving water up the column?

Day 4.5: Structures look different

Day 7: Waves visible in Omega (left)

Here is some initial timing information. This is not specific to this PR, but I thought I would post it here because the pressure gradient allows us to collect timing data for primitive equation tests for the first time. Using the branch hyungyukang/hkang/omega/new-test-branch with this PR fully merged in.

Summary of Omega vs MPAS-O timing with primitive equation terms.
Test is 1km baroclinic channel with 92k cells on 4 perlmutter nodes, 512 CPUs. Timing averaged over 10 RK4 timesteps. Timing includes the first time step for both models. The first time step is often 50% slower.

Timing without i/o:
Omega: 0.02658s/timestep (3% slower)
MPASO: 0.02579s/timestep

Full timing, with i/o, init, finalize for 10 RK4 timesteps:
Omega: 10.273s (2.3% slower)
MPASO: 10.041s

For reference, unsplit-RK4 dt=5s (left) versus split-AB2 dt=30s (right), both in MPAS-O, produce negligible differences. This is 1km resolution, day 4.

As @mark-petersen noted, I have also observed somewhat different behavior between Omega and MPAS-O in the baroclinic channel test at 4 km resolution.

Below are my results from the 4 km baroclinic channel test, along with a comparison between Omega and MPAS-Ocean. Waves appear in Omega by day 20, similar to what @mark-petersen reported at 1 km resolution.

The time series of the domain-integrated kinetic energy:

Omega and MPAS-O begin to diverge after about day 5, and Omega shows slightly higher total KE.

@mark-petersen and @hyungyukang, if you said, U missed it but are you using a second order pressure gradient in MPAS-O as well in your tests?

@mark-petersen and @hyungyukang, if you said, U missed it but are you using a second order pressure gradient in MPAS-O as well in your tests?

@xylar , I’ve been using pressure_and_zmid in MPAS-O, and I believe @mark-petersen has been doing the same (Please correct me if you’re using a different one.).

Okay, that's likely the closest equivalent. Sounds right.

I ran experiments with real bottom topography at QU240km resolution without forcing (kind of a spin-down configuration).

Omega blows up within 4 hours, producing NaN values in the variables, whereas MPAS-O remains stable through 20 days. When the PGF term is disabled, although this is not a physically appropriate configuration, Omega also remains stable through 20 days. The instability in Omega may be caused either by its interaction with other terms or by the PGF term itself, and this will need to be investigated further.

FYI, the domain integrated KE in MPAS-Ocean:

Merged locally into develop. Passes all cTests with gnu on perlmutter and Frontier, CPU and GPU on both.